Mapeando datos en R
Martín Alalú
29/04/2020
1 Intro datos geo
Hoy en día vemos mapas por todos lados. Si hay datos, entonces hay un mapa. De hecho se dice que el 80% de los datos tiene componente geo, o sea que casi todo puede ser ubicado en un mapa! Twitts, fotos de IG, pacientes, hospitales, comercios, noticias, etc, todo puede ser ubicado en un mapa.
El atributo que distingue a los datos georreferenciados es que representan ubicaciones exactas sobre la superficie de la Tierra. Es decir que cómo mínimo un dato georreferenciado tiene un par de coordenadas (X,Y) que representan la ubicación de un elemento en la superficie terrestre. Latitud y longitud es el sistema de coordenadas de referencia más usado hoy en día (CRS).
Otro concepto importante al momento de trabajar con datos geográficos es el de proyección cartográfica. Transformar una esfera (la tierra) en un plano (mapa) es todo un tema, no hay una única forma de hacerlo. La proyección cartográfica indica la serie de instrucciones que aplicamos para pasar a 2 dimensiones los distintos elementos de la esfera terrestre. Es por esto que hay muchísimas proyecciones cartográficas, y todas tienen ventajas y desventajas.
Al momento de manipular datos geográficos es sumamente importante saber con qué sistema de coordenadas de referencia y la proyección cartográfica fueron creados así podemos compararlos y trabajar con ellos adecuadamente.
En cuanto a tipos de datos geográficos encontramos 3:
- Puntos: Indica la posición exacta de un elemento (hospitales, reclamos, etc.). Está formado por sólo un par de coordenadas (X,Y)
- Líneas: Indica un recorrido (calles, ríos,etc.). Está formado por una secuencia de pares de coordenadas (X,Y) que no delimitan un área.
- Polígonos: Indica el área / superficie de un elemento (barrio, área hospitalaria, etc.). Está formado por una secuencia de pares de coordenadas que empiezan y terminan en el mismo punto, generando así un área cerrada, delimitada.
Para manipular datos geográficos en R vamos a usar la librería sf, una de las más populares en R y con muchísimas funciones disponibles. Al igual que con el resto de los paquetes de R lo instalamos usando install.packages(sf) y lo llamamos con library(sf).
library(sf)
library(tidyverse)2 Analisis espacial
El análisis espacial es un sub conjunto dentro del análisis de datos, que tiene como objetivo entender, describir, analizar, modelar, predecir la ubicación en el espacio de nuestra unidad de análisis. La premisa principal es que la ubicación per se es ya un elemento sumamente importante para comprender el fenómeno.
Conseguir datos geográficos es relativamente fácil hoy en día. Hay muchísimos fuentes de información que disponibilizan este tipo de datos:
En cuanto a formatos de archivos geográficos principalmente encontraremos 3: CSV, Shapefile (SHP) y GeoJson.
En esta ocasión trabajaremos con un archivo en formato GeoJson que tiene los radios censales de la Ciudad de Buenos Aires, descargado del portal de datos abiertos de la Ciudad. Los radios son unidades geográficas que agrupan en promedio 300 viviendas en las ciudades y se usan para relevar información en el terreno. En los censos los radios son la unidad geográfica mínima en la que se publica información referida a población, hogares y viviendas.
## Reading layer `caba_radios' from data source `https://github.com/martoalalu/clase-geo-salud/raw/master/data/caba_radios.geojson' using driver `GeoJSON'
## Simple feature collection with 3554 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -58.53092 ymin: -34.70574 xmax: -58.33455 ymax: -34.528
## Geodetic CRS: WGS 84Con solo cargar el archivo ya R nos da información valiosa sobre este archivo:
Simple feature collection with 3554 features and 8 fieldsya nos indica que el archivo en cuestión tiene 3554 observaciones y 8 columnas.geometry type: MULTIPOLYGONindica que el archivo geografico tiene multiples poligonos (en este caso radios censales)dimension: XYindica que tiene 2 dimensiones (o sea que no hay una tercera dimensión “Z”)bboxindica el bounding box, la caja que delimita y contiene a los datos en cuestión. Son los valores máximos y mínimos de las coordenadas.CRS: 4326indica el sistema de coordenadas de referencia en que están los datos. En este caso es el famoso WGS84.
Al fin y al cabo un dato geográfico es como cualquier dataframe solo que tiene un campo con la información de la geometría. Por eso vamos a ver inspeccionar qué tiene el dataframe radios.
head(radios)## Simple feature collection with 6 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -58.38593 ymin: -34.60846 xmax: -58.35054 ymax: -34.57865
## Geodetic CRS: WGS 84
## RADIO_ID BARRIO COMUNA POBLACION VIVIENDAS HOGARES HOGARES_NBI
## 1 1_1_1 RETIRO 1 336 82 65 19
## 2 1_12_1 SAN NICOLAS 1 341 365 116 25
## 3 1_12_10 SAN NICOLAS 1 296 629 101 1
## 4 1_12_11 SAN NICOLAS 1 528 375 136 7
## 5 1_12_2 SAN NICOLAS 1 229 445 129 16
## 6 1_12_3 SAN NICOLAS 1 723 744 314 104
## AREA_KM2 geometry
## 1 1.79899705 MULTIPOLYGON (((-58.37189 -...
## 2 0.01856469 MULTIPOLYGON (((-58.38593 -...
## 3 0.04438025 MULTIPOLYGON (((-58.37879 -...
## 4 0.36634000 MULTIPOLYGON (((-58.36733 -...
## 5 0.01836301 MULTIPOLYGON (((-58.38454 -...
## 6 0.03672540 MULTIPOLYGON (((-58.38154 -...summary(radios)## RADIO_ID BARRIO COMUNA POBLACION
## Length:3554 Length:3554 Length:3554 Min. : 0.0
## Class :character Class :character Class :character 1st Qu.: 646.2
## Mode :character Mode :character Mode :character Median : 786.0
## Mean : 813.2
## 3rd Qu.: 928.0
## Max. :3945.0
## VIVIENDAS HOGARES HOGARES_NBI AREA_KM2
## Min. : 0.0 Min. : 0.0 Min. : 0.00 Min. :0.004468
## 1st Qu.: 311.2 1st Qu.: 259.0 1st Qu.: 2.00 1st Qu.:0.018626
## Median : 377.0 Median : 310.0 Median : 6.00 Median :0.035548
## Mean : 401.4 Mean : 323.6 Mean : 19.35 Mean :0.057350
## 3rd Qu.: 462.0 3rd Qu.: 371.0 3rd Qu.: 23.00 3rd Qu.:0.062847
## Max. :1405.0 Max. :1093.0 Max. :403.00 Max. :3.804422
## geometry
## MULTIPOLYGON :3554
## epsg:4326 : 0
## +proj=long...: 0
##
##
## El dataset tiene un montón de información valiosa que podemos desplegar en un mapa. Vamos con el primero.
3 De la tabla al mapa
Bueno, ya basta de tablas, vamos al mapa! Vamos con una pregunta fácil. ¿Cómo se distribuye la población a lo largo y ancho de la Ciudad? ¿Donde hay más densidad poblacional?
Para hacerlo vamos a seguir utilzando la sintaxis de ggplot() pero indicando que ahora estamos trabajando con un dato geográfico, para eso vamos a usar geom_sf() como parámetro.
ggplot() +
geom_sf(data=radios)
Bien! Tal como podemos ver hay radios censales más grandes que otros, pero generalmente mantienen el mismo tamaño en este caso. Vamos con un primer mapa coloreando los radios censales según el barrio.
ggplot() +
geom_sf(data = radios, aes(fill = BARRIO))
Uuu! Tenemos muchísimas variables y no podemos ver el mapa! Saquemos la referencia.
ggplot() +
geom_sf(data = radios, aes(fill = BARRIO))+
guides(fill=FALSE)
Ahora bien, intentemos encontrar algunos patrones haciendo un mapa coroplético (choropleth) Vamos con un simple mapa que muestre la cantidad de población por radio así vemos las zonas más pobladas de la Ciudad. Para eso vamos a agregar nuevos parámetros a nuestro ggplot() y de paso agregamos un título.
ggplot() +
geom_sf(data = radios, aes(fill = POBLACION)) +
labs(title = "Población por radio censal",
subtitle = "Ciudad de Buenos Aires",
fill = "Habitantes")
Ahí va mejor. Pero las líneas grises de los límites de los radios son algo molestas, vamos a sacarlas indicandole que no queremos ningún color poniendo color = NA. Como estamos trabajando con mapas a veces la referencia de los ejes X e Y con la latitud y la longitud no nos interesa mucho, para sacarlas podemos directamente aplicar un theme que haga el mapa más “limpio”.
ggplot() +
geom_sf(data = radios, aes(fill = POBLACION), color = NA) +
labs(title = "Población por radio censal",
subtitle = "Ciudad de Buenos Aires",
fill = "Habitantes")+
theme_void()
Va quedando pero tampoco podemos sacar una conclusión clara sobre la distribución de la población en la Ciudad. Vemos algunos radios grandes que tienen mucha población (al sur en Villa Lugano o en el centro el de Chacarita) y otros que tienen poca (Bosques de Palermo y el Autódromo).
Vamos a aprovechar que una columna indica la superficie del radio en km2 y otra la población y vamos a crear una nueva columna que indique la densidad poblacional de cada radio, su población por km2. Seguimos aprovechando las bondades de tidyverse y usamos mutate.
Además vamos a cambiar la escala de cólores y usar [viridis](https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html), también muy usada, y que nos permitirá destacar mejor los patrones escondidos en los datos. El parámetro scale_fill_viridis_c() nos va a dejar hacer esto.
#Creamos nueva columna con densidad poblacional por barrio
radios <- radios %>%
mutate (densidad=POBLACION/AREA_KM2)
#Mapa según densidad de población usando viridis como escala de colores
ggplot() +
geom_sf(data = radios, aes(fill = densidad), color = NA) +
scale_fill_viridis_c() +
labs(title = "Densidad de población por radio censal",
subtitle = "Ciudad Autónoma de Buenos Aires",
fill = "hab/km2")+
theme_void()
Ahora sí podemos ver algunos patrones. El corredor Cabildo-Santa Fe, Rivadavia y Pueyrredón emergen como zonas más densamente pobladas al igual que varios barrios vulnerables (31, 1-11-14, 21, Soldati).
4 Coyuntura.. ¿¡o qué!? v2
Al igual que con cualquier dataframe podemos hacer agrupaciones de datos geográficos y “unir” polígonos según determinado criterio. Esto es sumamente útil cuando queremos pasar a trabajar con otras unidades geográficas.
Por ejemplo vamos a transformar este dataframe donde cada unidad es un radio censal a uno donde sean los barrios (o sea que pasaríamos de tener 3554 observaciones a 48). Para hacer esto usamos los conocidos group_by() y summarise().
barrios <- radios %>%
group_by(BARRIO) %>%
summarise(POBLACION = sum(POBLACION),
VIVIENDAS = sum(VIVIENDAS),
HOGARES = sum(HOGARES),
HOGARES_NBI = sum(HOGARES_NBI),
AREA_KM2 = sum(AREA_KM2))
head(barrios)## Simple feature collection with 6 features and 6 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -58.50294 ymin: -34.66295 xmax: -58.33455 ymax: -34.53199
## Geodetic CRS: WGS 84
## # A tibble: 6 x 7
## BARRIO POBLACION VIVIENDAS HOGARES HOGARES_NBI AREA_KM2
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AGRON… 13912 6262 5284 70 2.12
## 2 ALMAG… 131699 71216 58327 3404 4.05
## 3 BALVA… 138926 77981 60387 7122 4.34
## 4 BARRA… 89452 33058 31249 3850 7.96
## 5 BELGR… 126267 71363 54666 542 7.73
## 6 BOCA 45113 18133 16287 3460 5.02
## # … with 1 more variable: geometry <GEOMETRY [°]>Y veamos cómo queda el nuevo mapa de barrios.
ggplot() +
geom_sf(data = barrios, aes(fill = BARRIO), color = NA) +
labs(title = "Barrios de la Ciudad de Buenos Aires")+
theme_void()+
guides(fill=FALSE)
Siguiendo con la coyuntura vamos a usar datos de covid a nivel barrio en la Ciudad de Buenos Aires que publica el Gobierno de la Ciudad en su portal de datos abiertos BA Data.
Vamos a leer el archivo como hacemos con cualquier dataframe no-geográfico.
#Leemos el archivo
covid <- read.csv(file = "https://raw.githubusercontent.com/martoalalu/clase-geo-salud/clase-2021/data/covid_barrios_caba.csv", sep=",")
summary(covid)## X barrio comunitario contacto
## Min. : 1.00 Length:48 Min. : 1222 Min. : 111.0
## 1st Qu.:12.75 Class :character 1st Qu.: 3070 1st Qu.: 490.2
## Median :24.50 Mode :character Median : 5442 Median : 769.0
## Mean :24.50 Mean : 6798 Mean :1140.0
## 3rd Qu.:36.25 3rd Qu.: 7232 3rd Qu.:1198.0
## Max. :48.00 Max. :27124 Max. :4557.0
##
## en_investigacion importado trabajador_de_la_salud na
## Min. : 109.0 Min. : 1.000 Min. : 22.0 Min. : 1.00
## 1st Qu.: 389.5 1st Qu.: 2.000 1st Qu.: 126.2 1st Qu.: 6.00
## Median : 620.0 Median : 4.000 Median : 207.5 Median : 10.00
## Mean : 804.8 Mean : 6.979 Mean : 295.3 Mean : 19.94
## 3rd Qu.: 843.2 3rd Qu.: 8.000 3rd Qu.: 350.5 3rd Qu.: 19.50
## Max. :3184.0 Max. :46.000 Max. :1162.0 Max. :133.00
## NA's :1 NA's :1
## total
## Min. : 1513
## 1st Qu.: 4499
## Median : 6936
## Mean : 9065
## 3rd Qu.: 9618
## Max. :33743
## head(covid)## X barrio comunitario contacto en_investigacion importado
## 1 1 AGRONOMIA 1223 111 138 3
## 2 2 ALMAGRO 14296 2596 1830 17
## 3 3 BALVANERA 16852 3377 2799 20
## 4 4 BARRACAS 11660 3163 1449 2
## 5 5 BELGRANO 14555 1105 876 17
## 6 6 BOCA 6407 1175 733 1
## trabajador_de_la_salud na total
## 1 34 4 1513
## 2 928 41 19708
## 3 1042 66 24156
## 4 382 10 16666
## 5 453 62 17068
## 6 236 7 8559Vamos a intentar hacer un mapa que indique la cantidad de casos de covid por barrio. El problema es que propiamente dicho este dataframe no es geográfico ya que no tiene una columna con la geometría del barrio. Pero tenemos la columna localidad que indica el barrio por lo que podemos aprovechar y hacer un join con nuestro dataframe geográfico de barrios.
#Hacemos join
barrios <- left_join(barrios,covid,by = c("BARRIO" = "barrio"))
head(barrios)## Simple feature collection with 6 features and 14 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -58.50294 ymin: -34.66295 xmax: -58.33455 ymax: -34.53199
## Geodetic CRS: WGS 84
## # A tibble: 6 x 15
## BARRIO POBLACION VIVIENDAS HOGARES HOGARES_NBI AREA_KM2
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AGRON… 13912 6262 5284 70 2.12
## 2 ALMAG… 131699 71216 58327 3404 4.05
## 3 BALVA… 138926 77981 60387 7122 4.34
## 4 BARRA… 89452 33058 31249 3850 7.96
## 5 BELGR… 126267 71363 54666 542 7.73
## 6 BOCA 45113 18133 16287 3460 5.02
## # … with 9 more variables: geometry <GEOMETRY [°]>, X <int>, comunitario <int>,
## # contacto <int>, en_investigacion <int>, importado <int>,
## # trabajador_de_la_salud <int>, na <int>, total <int>Genial. Ahora tenemos un dataframe geográfico con cantidad de casos de covid por barrio. Vamos a hacer un nuevo mapa cloroplético con la tasa de covid cada 10.000 habitantes.
barrios <- barrios %>%
mutate(tasa_covid=(total/POBLACION)*10000)
ggplot(data=barrios) +
geom_sf(aes(fill = tasa_covid), color = NA) +
scale_fill_viridis_c() +
labs(title = "Casos covid por barrio cada 10000 habitantes",
subtitle = "Ciudad Autónoma de Buenos Aires",
fill = "Casos por cada 10000 habitantes")+
theme_void()
Flores emerge claramente como el barrio con la tasa más alta de covid por cada 10.000 habitantes. Núñez y Villa Real al oeste también. La zona sur parece la menos afectada por ahora.
Si queremos salir de la perspectiva geográfica podemos hacer un gráfico de barras como con cualquier dataframe.
ggplot(barrios) +
geom_bar(aes(x= reorder(BARRIO, tasa_covid), weight=tasa_covid)) +
coord_flip() 
5 Del polígono al punto
Por ahora estuvimos manejando un tipo de dato geográfico (polígonos) y haciendo un estilo de mapas (coropléticos), probemos usando datos que sean puntos y hagamos otros tipos de mapas.
Vamos a usar un dataset disponibilizado por el Ministerio de Salud del Gobierno Nacional que contiene la ubicación exacta todos los efectores de salud del país.
#Cargamos los datos
salud <- read.csv("https://github.com/martoalalu/clase-geo-salud/raw/master/data/refes-hospitales.csv", fileEncoding = 'UTF-8')
#Vemos la clase de objeto que es
class(salud)## [1] "data.frame"head(salud)## CODIGO GLN NOMBRE
## 1 5.038084e+13 NA PUESTO EL TORO (17)
## 2 5.038063e+13 NA PUESTO BARRO NEGRO (9) (EX LOTE ROSARIO DE RIO GRANDE)
## 3 5.026021e+13 NA KIN SPORT
## 4 5.126035e+13 NA LABORATORIO ANALISIS CLINICOS GEROSA
## 5 5.126035e+13 NA LABORATORIO ANALISIS CLINICOS
## 6 5.126035e+13 NA LABORATORIO ANATOMIA PATOLOGICO (SALVAGNO)
## CODIGO_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSIDT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## TIPOLOGIA
## 1 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 2 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 3 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 4 Establecimiento de salud sin internación de diagnóstico
## 5 Establecimiento de salud sin internación de diagnóstico
## 6 Establecimiento de salud sin internación de diagnóstico
## CODIGO_CATEGORIA_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSDIT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## CATEGORIA_TIPOLOGIA
## 1 Sin atención médica en forma periódica (menor a 3 veces por semana)
## 2 Con atención médica general por lo menos 3 días de la semana
## 3 Con atención médica diaria y con especialidades y/o otras profesiones
## 4 Laboratorio de Análisis Clínicos
## 5 Laboratorio de Análisis Clínicos
## 6 Laboratorio de Anatomía patológica
## CODIGO_DEPENDENCIA DEPENDENCIA ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA
## 1 21 Provincial Público 38
## 2 21 Provincial Público 38
## 3 23 Privado Privado 26
## 4 23 Privado Privado 26
## 5 23 Privado Privado 26
## 6 23 Privado Privado 26
## PROVINCIA CODIGO_DEPARTAMENTO DEPARTAMENTO CODIGO_LOCALIDAD
## 1 Jujuy 84 Susques 38084030000
## 2 Jujuy 63 San Pedro 38063090000
## 3 Chubut 21 Escalante 26021030018
## 4 Chubut 35 Futaleufú 26035030000
## 5 Chubut 35 Futaleufú 26035030000
## 6 Chubut 35 Futaleufú 26035030000
## LOCALIDAD CP DOMICILIO TE1
## 1 EL TORO 4411 AVENIDA LAVALLE ESQUINA 24 DE OCTUBRE
## 2 LA MENDIETA 4522 CARLOS GARDEL S/N° - LOTE BARRO NEGRO
## 3 COMODORO RIVADAVIA 3240 Avenida Rivadavia 633 0297-4475310
## 4 ESQUEL 9200 25 de Mayo 498 2945452050
## 5 ESQUEL 9200 Ameghino 1033 02945-452155
## 6 ESQUEL 9200 Belgrano 347
## TE2 TE3 TE4 LATITUD LONGITUD ESTADO_GEO
## 1 -23.20853 -66.82709 SINCONFIRMAR
## 2 -24.30814 -64.93175 SINCONFIRMAR
## 3 -45.86213 -67.48355 CONFIRMADO
## 4 294515503809 -42.91505 -71.31971 CONFIRMADO
## 5 -42.91577 -71.31803 CONFIRMADO
## 6 -42.91788 -71.31992 SINCONFIRMARsummary(salud)## CODIGO GLN NOMBRE CODIGO_TIPOLOGIA
## Min. :1.002e+13 Min. :7.798e+12 Length:28366 Length:28366
## 1st Qu.:5.006e+13 1st Qu.:9.990e+12 Class :character Class :character
## Median :5.050e+13 Median :9.991e+12 Mode :character Mode :character
## Mean :4.380e+13 Mean :9.868e+12
## 3rd Qu.:5.107e+13 3rd Qu.:9.992e+12
## Max. :5.394e+13 Max. :9.993e+12
## NA's :28241
## TIPOLOGIA CODIGO_CATEGORIA_TIPOLOGIA CATEGORIA_TIPOLOGIA
## Length:28366 Length:28366 Length:28366
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## CODIGO_DEPENDENCIA DEPENDENCIA ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA
## Min. :20.00 Length:28366 Length:28366 Min. : 2.00
## 1st Qu.:22.00 Class :character Class :character 1st Qu.: 6.00
## Median :23.00 Mode :character Mode :character Median :30.00
## Mean :22.52 Mean :38.81
## 3rd Qu.:23.00 3rd Qu.:66.00
## Max. :32.00 Max. :94.00
##
## PROVINCIA CODIGO_DEPARTAMENTO DEPARTAMENTO CODIGO_LOCALIDAD
## Length:28366 Min. : 1.0 Length:28366 Min. :6.028e+06
## Class :character 1st Qu.: 28.0 Class :character 1st Qu.:6.588e+09
## Mode :character Median : 70.0 Mode :character Median :3.004e+10
## Mean :157.5 Mean :3.749e+10
## 3rd Qu.:140.0 3rd Qu.:6.603e+10
## Max. :882.0 Max. :9.401e+10
##
## LOCALIDAD CP DOMICILIO TE1
## Length:28366 Length:28366 Length:28366 Length:28366
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## TE2 TE3 TE4 LATITUD
## Length:28366 Length:28366 Length:28366 Min. :-54.83
## Class :character Class :character Class :character 1st Qu.:-34.70
## Mode :character Mode :character Mode :character Median :-32.91
## Mean :-32.52
## 3rd Qu.:-28.47
## Max. :-21.94
## NA's :2240
## LONGITUD ESTADO_GEO
## Min. :-72.89 Length:28366
## 1st Qu.:-65.35 Class :character
## Median :-61.93 Mode :character
## Mean :-62.47
## 3rd Qu.:-58.67
## Max. :-53.65
## NA's :2240Lo cargamos como un dataframe normal no-geográfico pero al explorarlo vemos que tiene 2 columnas que sí precisan su ubicación exacta en la tierra, LATITUD y LONGITUD. Tenemos que hacerle entender a R que se trata de un dataframe geográfico, es fácil!
Primero filtramos y sólo nos quedamos con las observaciones que no son nulas ya que no podemos ubicar en un mapa objetos con latitud y longitud nulas. Luego con st_as_sf convertimos el dataframe a un objeto geográfico indicandole dónde están los datos de las coordenadas y le indicamos el sistema de referencia de coordenadas.
salud <- salud %>%
filter(!is.na(LONGITUD), !is.na(LATITUD)) %>%
st_as_sf(coords = c("LONGITUD", "LATITUD"), crs = 4326)
class(salud)## [1] "sf" "data.frame"head(salud)## Simple feature collection with 6 features and 23 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -71.31992 ymin: -45.86213 xmax: -64.93175 ymax: -23.20853
## Geodetic CRS: WGS 84
## CODIGO GLN NOMBRE
## 1 5.038084e+13 NA PUESTO EL TORO (17)
## 2 5.038063e+13 NA PUESTO BARRO NEGRO (9) (EX LOTE ROSARIO DE RIO GRANDE)
## 3 5.026021e+13 NA KIN SPORT
## 4 5.126035e+13 NA LABORATORIO ANALISIS CLINICOS GEROSA
## 5 5.126035e+13 NA LABORATORIO ANALISIS CLINICOS
## 6 5.126035e+13 NA LABORATORIO ANATOMIA PATOLOGICO (SALVAGNO)
## CODIGO_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSIDT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## TIPOLOGIA
## 1 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 2 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 3 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 4 Establecimiento de salud sin internación de diagnóstico
## 5 Establecimiento de salud sin internación de diagnóstico
## 6 Establecimiento de salud sin internación de diagnóstico
## CODIGO_CATEGORIA_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSDIT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## CATEGORIA_TIPOLOGIA
## 1 Sin atención médica en forma periódica (menor a 3 veces por semana)
## 2 Con atención médica general por lo menos 3 días de la semana
## 3 Con atención médica diaria y con especialidades y/o otras profesiones
## 4 Laboratorio de Análisis Clínicos
## 5 Laboratorio de Análisis Clínicos
## 6 Laboratorio de Anatomía patológica
## CODIGO_DEPENDENCIA DEPENDENCIA ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA
## 1 21 Provincial Público 38
## 2 21 Provincial Público 38
## 3 23 Privado Privado 26
## 4 23 Privado Privado 26
## 5 23 Privado Privado 26
## 6 23 Privado Privado 26
## PROVINCIA CODIGO_DEPARTAMENTO DEPARTAMENTO CODIGO_LOCALIDAD
## 1 Jujuy 84 Susques 38084030000
## 2 Jujuy 63 San Pedro 38063090000
## 3 Chubut 21 Escalante 26021030018
## 4 Chubut 35 Futaleufú 26035030000
## 5 Chubut 35 Futaleufú 26035030000
## 6 Chubut 35 Futaleufú 26035030000
## LOCALIDAD CP DOMICILIO TE1
## 1 EL TORO 4411 AVENIDA LAVALLE ESQUINA 24 DE OCTUBRE
## 2 LA MENDIETA 4522 CARLOS GARDEL S/N° - LOTE BARRO NEGRO
## 3 COMODORO RIVADAVIA 3240 Avenida Rivadavia 633 0297-4475310
## 4 ESQUEL 9200 25 de Mayo 498 2945452050
## 5 ESQUEL 9200 Ameghino 1033 02945-452155
## 6 ESQUEL 9200 Belgrano 347
## TE2 TE3 TE4 ESTADO_GEO geometry
## 1 SINCONFIRMAR POINT (-66.82709 -23.20853)
## 2 SINCONFIRMAR POINT (-64.93175 -24.30814)
## 3 CONFIRMADO POINT (-67.48355 -45.86213)
## 4 294515503809 CONFIRMADO POINT (-71.31971 -42.91505)
## 5 CONFIRMADO POINT (-71.31803 -42.91577)
## 6 SINCONFIRMAR POINT (-71.31992 -42.91788)Como ven el objeto salud pasó de ser un dataframe a un sf dataframe, o sea un dataframe geográfico. Y también podemos ver todas las características del objeto. El tipo de geometría es punto, tiene 2 dimensiones y el CRS es 4326 el mismo que el de radios.
Ahora podemos mapear. Y lo hacemos igual que con el dataset de barrios y radios.
ggplot()+
geom_sf(data=salud)
No tenemos los límites de Argentina pero claramente los puntos de los efectores parecen distribuirse en todo el país, es decir que no se limitan a aquellos que se ubican en Buenos Aires.
A no desesperarse! Una de las bondades de los datos geográficos es que podemos hacer consultas espaciales. En este caso vamos a pedirle a R que espacialmente nos filtre y seleccione todos los efectores de salud que caen dentro de algún barrio de Buenos Aires. Para eso usamos st_intersection una de las funciones espaciales de la librería sf.
#Nos quedamos solo con los efectores de salud que están dentro de barrios
salud_caba <- st_intersection(salud, barrios)La consulta espacial dio como resultado un nuevo dataframe que nos indica que en la Ciudad de Buenos Aires hay 1154 efectores de salud.
Ahora sí, volvamos al mapa. Al igual que en cualquier gráfico de ggplot() podemos combinar capas. Vamos a agregar una capa con los barrios detrás para tener la referencia.
ggplot()+
geom_sf(data=barrios)+
geom_sf(data=salud_caba)
A priori la distribución es parecida a la de densidad poblacional!
Ejercicio. Supongamos que querramos hacer un mapa que muestre la capacidad de respuesta del sistema de salud porteño, o sea sólo aquellos efectores que tengan internación. ¿Cómo lo hacemos? Va tip, el campo TIPOLOGÍA nos va a ayudar.
Volvamos al mapa, pero profundizando el análisis. Vamos a pintar cada efector según la columna ORIGEN_FINANCIAMIENTO el cual indica si es público o privado para ver cómo es la distribución espacial.
ggplot()+
geom_sf(data=barrios)+
geom_sf(data=salud_caba, aes(color=ORIGEN_FINANCIAMIENTO))
Para emprolijar podemos hacer facetar el mapa, al igual que con cualquier visualización de ggplot().
ggplot()+
geom_sf(data=barrios)+
geom_sf(data=salud_caba, aes(color=ORIGEN_FINANCIAMIENTO))+
facet_wrap(~ORIGEN_FINANCIAMIENTO)
Vamos con otro ejercicio ¿Cómo sería si quisieramos hacer un nuevo mapa choropleth con la cantidad de efectores de salud con internación por barrio?
6 Más allá del choropleth
Hagamos otros tipos de mapas. Tener ubicaciones exactas, es decir objetos que son puntos, nos permiten hacer otros tipos de mapas como de calor (heatmap), de burbujas (bubblemap) o de densidad (densitymap), entre otros.
Dejemos de usar un rato la librería sf y vamos a considerar nuestro dataframe de efectores de salud en Buenos Aires como cualquier otro no geográfico.
Antes con st_as_coord() habíamos transformado las columnas de latitud y longitud en una geometry ahora vamos a reestablecer esas columnas.
#Reestablecemos las columnas lat y long con sus coordenadas
salud_caba$long <- st_coordinates(salud_caba)[,1]
salud_caba$lat <- st_coordinates(salud_caba)[,2]
#Le decimos a R que salud_caba_intern deja de ser un objeto geográfico
st_geometry(salud_caba) <- NULL
class(salud_caba)## [1] "data.frame"Ahora salud_caba es un dataframe normal que tiene una columna latitud y otra longitud que expresan la posición de los elementos en la superficie terrestre. O sea que si hacemos un
Volvamos a ggplot() y hagamos un geom_point() tradicional.
ggplot() +
geom_point(data = salud_caba, aes(x = long, y = lat,color = ORIGEN_FINANCIAMIENTO))
Ven que queda igual que si fuera un objeto geográfico! Aprovechemos nuevas funciones de ggplot() para hacer nuevos gráficos que “simulen” ser un mapa. Un mapa de densidad de puntos que muestre las zonas con mayor concentración de efectores.
ggplot() +
geom_bin2d(data = salud_caba, aes(x = long, y = lat), bins = 50) +
scale_fill_viridis_c()+
labs(title = "Concentración de efectores de salud",
fill = "Cantidad")
A todo esto siempre estuvimos trabajando con puntos y polígonos asumiendo que el fondo es Buenos Aires. Pero podemos agregar una capa base que tenga el mapa de Buenos Aires.
Para eso vamos a usar las librerías ggmap y osmdata, las cuales instalamos con el tradicional install.packages("NOMBRE_LIBRERIA") y agregamos con library(NOMBRE_LIBRERIA). Con estas nuevas librería podemos agregar capas de fondo muy fácilmente.
Básicamente lo primero que tenemos que hacer es indicarle a R cuál es el limite, el boundin gbox del mapa que queremos descargarnos. Usamos la función getbb y le decimos que es “Ciudad Autónoma de Buenos Aires” ya que es el nombre oficial. Si estuvieramos haciendo un mapa de por ejemplo Rosario le indicaríamos “Rosario, Santa Fe, Argentina”.
Luego, con get_stamenmap le indicamos especialmente que nos descargue el mapa, en este caso el tipo de mapa que nos estamos descargando es el toner-background, pero si quieren otro tipo pueden chequear la página de Stamen Maps.
Luego, en vez de usar ggplot() vamos a usar ggmap() pero no se preocupen la sintaxis sigue siendo la misma!
library(ggmap)
library(osmdata)
#Indicamos cual es la ciudad de fondo que queremos
bbox <- getbb("Ciudad Autónoma de Buenos Aires,Argentina")
#Descargamos el mapa
CABA <- get_stamenmap(bbox = bbox,
maptype = "terrain",zoom=12)
# Mapa solo
ggmap(CABA)
Nos faltan los datos pero ahora sabemos con certeza que el fondo corresponde a Buenos Aires.
Si no nos gusta este mapa podemos probar con otro, por ejemplo con toner-background.
#Descargamos el mapa
CABA <- get_stamenmap(bbox = bbox,
maptype = "toner-background",zoom=12)
ggmap(CABA)
Ahora si, agregemos los datos!
ggmap(CABA)+
geom_bin2d(data = salud_caba, aes(x = long, y = lat), bins = 50) +
scale_fill_viridis_c()
Excelente! Ahora probemos con un mapa que muestre la misma información pero simulando curvas de nivel.
ggmap(CABA) +
stat_density2d(data = salud_caba, aes(x = long, y = lat, fill=stat(level)),geom="polygon")+
scale_fill_viridis_c()